
hdfFilePathy = img;

expp = ExpParams();
cons = Constants();
img_time = 100e-6; % Length of imaging pulse (100 us)

% y-signal (in camera frame, not using expt axes)
yInt = sum(y,2);
xInt = sum(y,1);

yC = 1:1040;
xC = 1:1392;

yMax = find(yInt == max(yInt),1);
xMax = find(xInt == max(xInt),1);

ampGuessY = max(yInt) - min(yInt);
ampGuessX = max(xInt) - min(xInt);

options = optimset('Display','off','TolFun',1e-16,'TolX',1e-16);
pY = fminsearch(@(p)FitGauss(p,yC,yInt),[ampGuessY yMax 1500 min(yInt)],options);
pX = fminsearch(@(p)FitGauss(p,xC',xInt),[ampGuessX xMax 1500 min(xInt)],options);

% Estimate number of atoms using camera
amp = (pX(1) + pY(1))/2; % Take mean amplitude of two Gaussian fits. 
capture_fraction = 0.5*(1 - sqrt(1 - (expp.na_telec)^2)); % fraction of total light captured by lens ( ~ 0.5 %)
total_num_electrons = amp*sqrt(2*pi)*mean([pX(3) pY(3)]);
total_photons = total_num_electrons/(expp.qe_std*expp.trans780*expp.trans_telec*capture_fraction);
num_atoms = round(2*total_photons/(img_time*cons.Gamma)); % num_atoms = (num_photons/time) / (scattering rate)

figure(1)
subplot(3,3,[1 2 4 5])
imagesc(y);
axis equal
xlims = xlim;
ylims = ylim;

subplot(3,3,[3 6])
plot(yInt,yC,gauss(pY,yC),yC);
set(gca,'YDir','reverse');
ylim(ylims);

subplot(3,3,[7 8])
plot(xC,xInt,xC,gauss(pX,xC));
xlim(xlims);

% Give the width as the std deviation in um
px_size = 6.45; % um
mag = 1/2; % magnification
pY(3) = (pY(3))*px_size/mag;
pX(3) = (pX(3))*px_size/mag;

h = subplot(3,3,9,'Visible', 'off');
cla(h);
text(0, .5, sprintf('%s\n%s\n%s\n%s', ...
        ['v-width (um): ' num2str(pX(3))], ['h-width (um): ' num2str(pY(3))], ...
        ['peak(v,h) = (' num2str(round(pX(2))) ',' num2str(round(pY(2))) ')'], ...
        ['Na = ' num2str(num_atoms,'%e')]), ...
        'VerticalAlignment', 'middle', ...
        'FontName', 'Courier New', ...
        'FontWeight', 'bold', ...
        'FontSize', 20);
